#version 120
#extension GL_EXT_gpu_shader4 : enable //access to integer and bit operation

const float voxelMax = 255.0;

uniform sampler2D backfacetexture;
uniform sampler3D volumetexture;

uniform usampler3D nodePool3D;
uniform sampler3D blockPool;


uniform sampler1D tranf;

uniform float stepsize;
uniform float begin;
uniform float isoValue;

varying vec4 vposition;
varying vec3 eposition;
uniform vec4 lightposition;


const float shininess = 5;
const float diffuse = 0.0;
const float specular = 1.0;

/////////////////////////////////////////////////////////

vec4 DirectVolumeRendering(vec3 front, vec3 back, float upperThreshold, float threshold, float offset)
{
	vec3 dir;
	dir = back - front;
	float len = length(dir);
	vec3 norm_dir = normalize(dir);
	//float iso_value = 0.1;

	vec3 rayIncrement = stepsize*norm_dir;

	vec3 pos = front;
	vec4 value=vec4(0,0,0,0);
	// accumlate value
	vec4 result = vec4(0,0,0,0);//RGBA
	float raylength = 0;

    float scale = 1.0/(upperThreshold - threshold);
    float res = 0;
    
	while(result.a < 1.0 && raylength < len)
	{
		float factor = 0.5;

	
		value = texture3D(blockPool, front+raylength*norm_dir);
		raylength += stepsize;

		vec4 texValue = texture1D(tranf, value.a);

		result.rgb += (1.0-result.a)*texValue.a*texValue.rgb;
		result.a += (1.0-result.a)*texValue.a;

	}
	return result;

}




//====node related======================================================

// GLSL support unsigned int (uint), not sure if uint can be compiled on ATI cards
struct nodeTexel
{
	uint first;
	uint second;
};
//###############Please Update this######################
//const !!!!!! fuck!!!!!! I need to change this when pool changes 
//head - 128
//uint nodePoolDim = 2u;
//uint blockPoolDim = 2u;
//uint blockSize = 129u;

//head - 64
//const uint nodePoolDim = 6u;
//const uint blockPoolDim = 4u;
//const float invBlockPoolDim = 1.0/float(blockPoolDim);
//const uint blockSize = 65u;
//const int octreeDepth = 2;
//head - 32
//const uint nodePoolDim = 10u;
//const uint blockPoolDim = 8u;
//const float invBlockPoolDim = 1.0/float(blockPoolDim);
//const uint blockSize = 33u;
//const int octreeDepth = 3;
//const uint nodePoolDim = 8u;
//const uint blockPoolDim = 6u;
//const float invBlockPoolDim = 1.0/float(blockPoolDim);
//const uint blockSize = 33u;
//const int octreeDepth = 3;
//beetle 64
//const uint nodePoolDim = 10u;
//const uint blockPoolDim = 7u;
//const float invBlockPoolDim = 1.0/float(blockPoolDim);
//const uint blockSize = 65u;
//const int octreeDepth = 4;

//beetle 32
//const uint nodePoolDim = 14u;
//const uint blockPoolDim = 11u;
//const float invBlockPoolDim = 1.0/float(blockPoolDim);
//const uint blockSize = 33u;

// 1G
const uint nodePoolDim = 28u;
const uint blockPoolDim = 23u;
const float invBlockPoolDim = 1.0/float(blockPoolDim);
const uint blockSize = 33u;
//const uint nodePoolDim = 50u;
//const uint blockPoolDim = 43u;
//const float invBlockPoolDim = 1.0/float(blockPoolDim);
//const uint blockSize = 17u;
//####################################################

//Utility
nodeTexel intTexture3D(usampler3D tex, uint z, uint y, uint x, uint maxZ, uint maxY, uint maxX)
{
	//maxZ (1:maxZ) number of element = maxZ
	// z = (0: max-1)
	float indexX = (float(x)+0.5/float(maxX))/float(maxX);
	float indexY = (float(y)+0.5/float(maxY))/float(maxY);
	float indexZ = (float(z)+0.5/float(maxZ))/float(maxZ);//TODO: convert to a vec operation

	uvec4 uSample = texture3D(tex, vec3(indexX, indexY, indexZ));

	nodeTexel node;
	node.first = uSample.r;
	node.second = uSample.a;
	return node;
}


uvec3 getNodeClusterPointer( nodeTexel node)
{
	uvec3 cPointer=uvec3(0u,0u,0u);
	cPointer[0] = node.first & uint(0x3FF) ;
	cPointer[1] = (node.first>>10u) & uint(0x3FF);
	cPointer[2] = (node.first>>20u) & uint(0x3FF);
	return cPointer;
}
nodeTexel getChildrenNode(uvec3 clusterPointer, uint k, uint j, uint i)
{
	nodeTexel node = intTexture3D(nodePool3D, clusterPointer[2]*2u+k, clusterPointer[1]*2u+j, clusterPointer[0]*2u+i,
													nodePoolDim,nodePoolDim,nodePoolDim);

	return node;
}

uvec3 getBlockPointer( nodeTexel node) 
{
	uvec3 bPointer=uvec3(0u,0u,0u);
	bPointer[0] = node.second & uint(0x3FF) ;
	bPointer[1] = uint( (node.second >> 10u) & uint(0x3FF));
	bPointer[2] = uint( (node.second >> 20u) & uint(0x3FF));
	return bPointer;
}

bool hasBlock( nodeTexel node)
{
	bool block = ((node.first>>30u &  uint(0x1)) >  uint(0))? true:false;
	return block;
}
bool isSubDividable( nodeTexel node)
{
	bool dividable =  ((node.first>>31u &  uint(0x1)) >  uint(0) )? true:false;
	return dividable;
}

//quick test
//function checked
//nodeTexel getChildrenNode(uvec3 clusterPointer, uint k, uint j, uint i) //##
//bool isSubDividable( nodeTexel node) //##
//bool hasBlock( nodeTexel node) //##
//uvec3 getNodeClusterPointer( nodeTexel node) //##
//uvec3 getBlockPointer( nodeTexel node) //##
vec3 OctreeCoord(vec3 pos)
{ 
	uvec3 nodeClusterPointer=uvec3(0u,0u,0u);
	uvec3 childrenIndex; //@@@@you can't write childrenIndex = uvec3(2*pos) !!!!!!!!!!WTF@@@@

	nodeTexel currentNode;
	currentNode.first = 1u<<31u; currentNode.second = 0u;

	uint octreeLevel = 0u;
	//float scaleFactor = 1.0; //@@ the concept of scale is wrong @@!!!!!!
	
	bool SubD=true;
	while(SubD) //@@can't write isSubDividable(currentNode);@@!!!WTF
	{
		//decending
		childrenIndex.x = uint((2*pos.x));
		childrenIndex.y = uint((2*pos.y));
		childrenIndex.z = uint((2*pos.z));
		pos = fract(2*pos);
		
		currentNode=getChildrenNode( getNodeClusterPointer(currentNode), childrenIndex.z, childrenIndex.y, childrenIndex.x);
		SubD = isSubDividable( currentNode );

		//scaleFactor *= 0.5;
	}

	//return vec3(childrenIndex);

	//calculate blockoffset in the pool
	uvec3 blockPointer=getBlockPointer(currentNode);

	vec3 offset = (invBlockPoolDim)*vec3(blockPointer);
	if(!hasBlock(currentNode))
		return vec3(0.0,0.0,0.0);

	//boundary shrink factor
	float shrinkFactor = (float(blockSize)-1)/(float(blockSize));
	vec3 shrinkOffset = vec3(0.5/float(blockSize))*invBlockPoolDim;
	//scale the pos
	return offset + (pos*invBlockPoolDim*shrinkFactor + shrinkOffset);
	//return offset + pos*invBlockPoolDim;
}

////////////////////////
struct range
{
	float tNear;
	float tFar;
};

struct BlockBox
{
	vec3 boxMin;
	vec3 boxMax;
	bool hasBlock;
	uvec3 blockPointer;
	float sampleRate;
	float singleColor;
};


range intersectBox(vec3 entryPoint , vec3 rayDir, vec3 boxmin, vec3 boxmax)//, float *tnear, float *tfar)
{
    // compute intersection of ray with all six bbox planes
    vec3 invR = vec3(1.0) / rayDir;
    vec3 tbot = invR * (boxmin - entryPoint);
    vec3 ttop = invR * (boxmax - entryPoint);

    // re-order intersections to find smallest and largest on each axis
    vec3 tmin = min(ttop, tbot);
    vec3 tmax = max(ttop, tbot);

    // find the largest tmin and the smallest tmax
    float largest_tmin = max(max(tmin.x, tmin.y), tmin.z);
    float smallest_tmax = min(min(tmax.x, tmax.y), tmax.z);

	range tMinMax;
	tMinMax.tNear = largest_tmin;
	tMinMax.tFar = smallest_tmax;

	return tMinMax;
}

//vec4 accumulate(vec3 rayO, vec3 rayDir, vec4 color, float tStart, float tEnd, float sampleRate)
//{
//	float  t=tStart;
//	vec4 value, texValue;
//		
//	while(color.a < 1.0 && t < tEnd)
//	{
//		value = texture3D(blockPool, rayO + t*rayDir);
//		texValue = texture1D(tranf, value.a);
//
//		color.rgb += (1.0-color.a)*texValue.a*texValue.rgb;
//		color.a += (1.0-color.a)*texValue.a;
//		t += sampleRate;
//
//	}
//
//	return color;
//}



BlockBox findCurrentBlockAndBoundingBox(vec3 pos)
{
	vec3 boxmin = vec3(0,0,0);
	vec3 boxmax;// = vec3(1,1,1);
	float boxsize = 1.0;//vec3(1,1,1);

	uvec3 nodeClusterPointer=uvec3(0u,0u,0u);
	uvec3 childrenIndex; //@@@@you can't write childrenIndex = uvec3(2*pos) !!!!!!!!!!WTF@@@@

	nodeTexel currentNode;
	currentNode.first = 1u<<31u; currentNode.second = 0u;

	
	float scaleFactor = 1.0;
	
	bool SubD=true;
	while(SubD) //@@can't write isSubDividable(currentNode);@@!!!WTF
	{
		//decending

		childrenIndex.x = uint((2*pos.x));
		childrenIndex.y = uint((2*pos.y));
		childrenIndex.z = uint((2*pos.z));

		boxsize *= 0.5;
		boxmin += vec3(childrenIndex)*boxsize;
		//scaleFactor *= 0.5;

		pos = fract(2*pos);
		
		currentNode=getChildrenNode( getNodeClusterPointer(currentNode), childrenIndex.z, childrenIndex.y, childrenIndex.x);
		SubD = isSubDividable( currentNode );

		//scaleFactor *= 0.5;
	}
	uvec3 blockPointer=getBlockPointer(currentNode);
	boxmax = boxmin+vec3(boxsize);

	BlockBox renderBlock;
	renderBlock.singleColor = float(currentNode.second)/voxelMax;
	renderBlock.hasBlock=hasBlock(currentNode);
	renderBlock.boxMin = boxmin;
	renderBlock.boxMax = boxmax;
	renderBlock.blockPointer = blockPointer;
	
	return renderBlock;
}
	////calculate blockoffset in the pool
	//uvec3 blockPointer=getBlockPointer(currentNode);

	//vec3 offset = (invBlockPoolDim)*vec3(blockPointer);
	//if(!hasBlock(currentNode))
	//	return vec3(0.0,0.0,0.0);

	////boundary shrink factor
	//float shrinkFactor = (float(blockSize)-1)/(float(blockSize));
	//vec3 shrinkOffset = vec3(0.5/float(blockSize))*invBlockPoolDim;
	////scale the pos
	//return offset + (pos*invBlockPoolDim*shrinkFactor + shrinkOffset);

//vec4 accumulate(vec3 rayO, vec3 rayDir, vec4 color, float tStart, float tEnd, float sampleRate)

//ray and dir is alway object coodinate
vec4 accumulateSingleColorBlock(vec4 color,BlockBox renderBlock, range r)
{
	//float invboxsize =1.0/( renderBlock.boxMax.r - renderBlock.boxMin.r);
	//float shrinkFactor = (float(blockSize)-1)/(float(blockSize))*invboxsize*invBlockPoolDim;
	////@@Need 
	//vec3 shrinkOffset = vec3(0.5/float(blockSize))*invBlockPoolDim;
	float tStart = r.tNear;
	float tEnd =  tStart+(r.tFar-r.tNear);
	//accumlate
	float  t=tStart;
	vec4 value, texValue;
	vec3 pos;	
	texValue = texture1D(tranf, renderBlock.singleColor);

	//TODO analytical integrate
	//safegard
	int i=0;
	while(color.a < 1.0 && t < tEnd)
	{
		if(i>200)
			break;
		i++;
		//pos = blockO+t*rayDir;
		//value = texture3D(blockPool, offset + pos*shrinkFactor + shrinkOffset);
		
		color.rgb += (1.0-color.a)*texValue.a*texValue.rgb;
		color.a += (1.0-color.a)*texValue.a;
		t += stepsize;

	}

	color.rgb += (1.0-color.a)*texValue.a*texValue.rgb;
	color.a += (1.0-color.a)*texValue.a;

	return color;

}

vec4 accumulateBlock(vec3 rayO, vec3 rayDir, vec4 color,BlockBox renderBlock, range r)
{
	//warping the coordinate
	uvec3 blockPointer = renderBlock.blockPointer;
	vec3 blockO = rayO - renderBlock.boxMin;
	vec3 offset = (invBlockPoolDim)*vec3(blockPointer);

	float invboxsize =1.0/( renderBlock.boxMax.r - renderBlock.boxMin.r);
	float shrinkFactor = (float(blockSize)-1)/(float(blockSize));
//	float shrinkFactor = (float(blockSize)-1)/(float(blockSize))*invboxsize*invBlockPoolDim;
	//@@Need 
	vec3 shrinkOffset = vec3(0.5/float(blockSize))*invBlockPoolDim;
	
	//vec3 invR = vec3(1.0) / rayDir;
	//float tStart = r.tNear+(shrinkOffset*rayDir);
	float tStart = r.tNear;//r.tNear;
	//float tEnd = tStart+shrinkFactor*(r.tFar-r.tNear);
	float tEnd =  r.tFar;
	if(r.tNear>r.tFar)
		return color;

	//accumlate
	float  t=0;
	//float t+=residual;
	vec4 value, texValue;
	vec3 pos;	

	//safeguard
	int i=0;
	while(color.a < 1.0 && t < tEnd)
	{
		i++;
		if(i>200)
			break;

		pos = blockO+t*rayDir;
		//return offset + (pos*invBlockPoolDim*shrinkFactor + shrinkOffset);
		value = texture3D(blockPool, offset + pos*shrinkFactor*invboxsize*invBlockPoolDim + shrinkOffset);
		texValue = texture1D(tranf, value.a);

		color.rgb += (1.0-color.a)*texValue.a*texValue.rgb;
		color.a += (1.0-color.a)*texValue.a;
		t += stepsize*shrinkFactor;

	}

	//visualize block Pointer
	//color += vec4(offset,1)*0.2;

	return color;
}

vec4  octreeRestartRayCasting(vec3 start, vec3 backPosition)
{
	vec4 color=vec4(0.0,0.0,0.0,0.0);
	vec3 rayDir = normalize(backPosition - start);
	float tMax = length(backPosition - start);
	float t = 0;
	
	int i=0;
	vec3 rayO;
	BlockBox renderBlock;
	range r ;
	while(t<tMax-0.001 && color.a <1.0)
	//for(int i=0; i<3; i++)
	{
		//ray decending, find the bounding box for current block
		rayO = start + t*rayDir;
		renderBlock = findCurrentBlockAndBoundingBox(rayO+0.00001*rayDir);
		//calculate ray box intersection

		r = intersectBox(rayO , rayDir, renderBlock.boxMin, renderBlock.boxMax);

		if(i==40)
		{
			//color = vec4(1,0,0,1);
			break;
		}
		i++;

		//accumulate volume
		if(renderBlock.hasBlock)
		{
			//color = vec4(rayO+r.tNear*rayDir,1);
			color = accumulateBlock(rayO, rayDir, color,renderBlock, r);
			//color = vec4(1,0,0,1);
		}
		else
		{
			color = accumulateSingleColorBlock(color,renderBlock,r);
		}

		//advance ray
		t += (r.tFar-r.tNear); 
	}

	return color;
}

vec4 DirectVolumeRenderingOctree(vec3 front, vec3 back, float upperThreshold, float threshold, float offset)
{
	vec3 dir;
	dir = back - front;
	float len = length(dir);
	vec3 norm_dir = normalize(dir);
	//float iso_value = 0.1;

	vec3 rayIncrement = stepsize*norm_dir;

	vec3 pos = front;
	vec4 value=vec4(0,0,0,0);
	// accumlate value
	vec4 result = vec4(0,0,0,0);//RGBA
	float raylength = 0;
//####
	//return vec4(OctreeCoord(pos+norm_dir*0.1),1.0);
//####

    float scale = 1.0/(upperThreshold - threshold);
    float res = 0;
    
	for(int i=0; i<3000; i++)
	{

		//value = texture3D(blockPool, pos);
		//if(pos.y>isoValue*1.6)
			value = texture3D(blockPool, OctreeCoord(pos));
		//else
		//	value = vec4(0,0,0,0);
	
		pos += rayIncrement;
		raylength += length(rayIncrement);

			vec4 texValue = texture1D(tranf, value.a);

			result.rgb += (1.0-result.a)*texValue.a*texValue.rgb;
			result.a += (1.0-result.a)*texValue.a;

			if (result.a >= 0.9 || raylength > len ) 
			{
				break;
			}
	}
	return result;

}

void main()
{
	vec2 texposition = ((vposition.xy/vposition.w) + 1.0)/2.0;

	vec3 start = gl_TexCoord[0].xyz;
	vec3 backposition = texture2D(backfacetexture, texposition).rgb;

	float upperThreshold = 1.0;
	float threshold = 0.1;
	vec4 result;
	if(isoValue<0.5)
		result = octreeRestartRayCasting(start, backposition);
		//result = DirectVolumeRenderingOctree(start, backposition, upperThreshold, threshold, 0);
	else
		result = octreeRestartRayCasting(start, backposition);
		//result = DirectVolumeRendering(start, backposition, upperThreshold, threshold, 0);
	gl_FragColor = (1-result.a)*vec4(0.2, 0.2, 0.2, 1) + result.a * result;

	//uvec3 p = uvec3(0u,0u,1u);
	//uvec3 color = getBlockPointer(getChildrenNode(p,1u,1u,1u));
	//gl_FragColor = vec4(vec3(color)/3.0,1);
	//gl_FragColor = vec4(float(isSubDividable(getChildrenNode(p,0u,0u,1u))),0,0,1);
}

//function checked
//nodeTexel getChildrenNode(uvec3 clusterPointer, uint k, uint j, uint i) //##
//bool isSubDividable( nodeTexel node) //##
//bool hasBlock( nodeTexel node) //##
//uvec3 getNodeClusterPointer( nodeTexel node) //##
//uvec3 getBlockPointer( nodeTexel node) //##
